EQUATION SOLVING 
Cross reference to related applications 

This application is a continuation in part of International Application 
5 PCT/GB03/001568, filed April 10, 2003, which claims priority to Great Britain 
Application No. GB 0208329.3, filed April 11, 2002, the contents of each of which 
are incorporated herein by reference. 

Field of Invention 

10 The present invention relates to systems and methods for solving systems of linear 
equations. 

Background of Invention 

Systems of linear equations occur frequently in many branches of science and 
15 engineering. Effective methods are needed for solving such equations. It is desirable 
that systems of linear equations are solved as quickly as possible. 

A system of linear equations typically comprises N equations in N unknown 
variables. For example, where N = 3 an example system of equations is set out 
20 below: 

15x + 5.y-2z = 15 
5x + 1 ly + 4z = 47 
-2x + 4j; + 9z = 51 

In this case, it is necessary to find values of jc, y, and z which satisfy all three 
25 equations. Many methods exist for finding such values of x, y and z and in this case 
it can be seen that x = 1, y = 2 and z = 5 is the unique solution to the system of 
equations. 

In general terms, existing methods for solving linear equations can be categorised as 
30 either direct methods, or iterative methods. Direct methods attempt to produce an 

1 



exact solution by using a finite number of operations. Direct methods however 
suffer from a problem in that the number of operations required is often large which 
makes the method slow. Furthermore, some implementations of such methods are 
sensitive to truncation errors. 

5 

Iterative methods solve a system of linear equations by repeated refinements of an 
initial approximation until a result is obtained which is acceptably close to the 
accurate solution. Each iteration is based on the result of the previous iteration, and, 
at least in theory, each iteration should improve the previous approximation. 
10 Generally, iterative methods produce an approximate solution of the desired accuracy 
by yielding a sequence of solutions which converge to the exact solution as the 
number of iterations tends to infinity. 

It will be appreciated, that systems of linear equations are often solved using a 
15 computer apparatus. Often, equations will be solved by executing appropriate 
program code on a microprocessor. In general terms, a microprocessor can represent 
decimal numbers in fixed point or floating point form. 

In fixed point form, a binary number of P bits comprises a first part of length q and a 
20 second part of length r. The first part represents the whole number part, and the 
second part represents the fractional part. In general terms, arithmetic operations can 
be relatively efficiently implemented for fixed point numbers. However, fixed point 
numbers suffer from problems of flexibility given that the position of the decimal 
point is fixed and therefore the range of numbers which can be accurately 
25 represented is relatively small given that overflow and round off errors regularly 
occur. 

Floating point numbers again in general terms comprise two parts. A first part, 
known as the exponent, and a second part known as the mantissa. The mantissa 
30 represents the binary number, and the exponent is used to determine the position of 
the decimal point within that number. 
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For example considering an eight bit value 00101011 where the first bit represents 
sign, the subsequent three bits represent the exponent, and the final four bits 
represent the mantissa, this value is analysed as follows. The first bit represents sign, 
and is interpreted such that the number is positive if the sign bit is equal to '0' and 
5 negative if the sign bit is equal to '1\ In this case, given that the first bit is '0', the 
number is determined to be positive. The mantissa (1011) is written out and a 
decimal point is placed to its left side as follows: 

.1011 

10 

The exponent is then interpreted as an integer, that is 010 is interpreted as the integer 
value 2 (given that the first bit of the exponent field is a sign bit which determines 
the direction in which the decimal point should move). In this case, the decimal 
point is moved to the right two bits to give: 

15 

10.11 

Which represents a value of two and three quarters. 

20 Although floating point numbers give considerable benefits in terms of their 
flexibility, arithmetic operations involving floating point numbers are inherently 
slower than corresponding operations on fixed point numbers. Therefore, where 
possible, the benefits of speed associated with fixed point numbers should be 
exploited. 

25 

When considering the implementation of an algorithm in hardware, its efficiency is a 
prime concern. Many algorithms for the solution of linear equations involve 
computationally expensive division and/or multiplication operations. These 
operations should, where possible be avoided, although this is often not possible with 
30 known methods for solving linear equations. 
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Many systems of linear equations have sparse solutions. In such cases the number of 
iterations required to solve the system of equations should be relatively low. 
However, this does not occur with some known methods. 



5 In summary there is a need for a more efficient algorithm which can be used to solve 
systems of linear equations. 

It is an object of the present invention to obviate or mitigate at least one of the 
problems identified above. 

10 

Summary of the Invention 

According to a first aspect of the present invention, there is provided a method for 
solving a system of N linear equations in N unknown variables, the method 
comprising: 

1 5 (a) storing an estimate value for each unknown variable; 

(b) initialising each estimate value to a predetermined value; 

(c) for each estimate value: 

(i) determining whether a respective predetermined condition 
is satisfied; and 

20 (ii) updating the estimate if and only if the respective 

predetermined condition is satisfied; and 

(d) repeating step (c) until each estimate value is sufficiently close to an 
accurate value of the respective unknown variable. 



25 Thus, the invention provides a method for solving linear equations in which 
estimates for solutions of the equations are updated only if a predetermined condition 
is satisfied. The predetermined condition is preferably related to convergence of the 
method. Therefore such an approach offers considerable benefits in terms of 
efficiency, given that updates are only carried out when such updates are likely to 

30 accelerate convergence. 
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According to a second aspect of the present invention, there is provided a method for 
solving a system of N linear equations in N unknown variables, the method 
comprising: 

(a) storing an estimate value for each unknown variable; 

(b) initialising each estimate value to a predetermined value; 

(c) attempting to update each estimate value using a scalar value d\ 

(d) updating the scalar value if no updates are made in step (c); and 

(e) repeating step (c) and step (d) until each estimate value is sufficiently 
close to an accurate value of the respective unknown variable. 

By updating the scalar value in accordance with the second aspect of the present 
invention, it has been discovered that benefits of efficiency are obtained. 



According to a third aspect of the present invention, there is provided a method for 
1 5 solving a system of N linear equations in N unknown variables of the form: 



Rh = /? 



the method comprising: 
20 generating a quadratic function of the form: 

J W) = ZZ R (^ w ) h ( w ) h ( w )- 2 Z^( w ) h (") ; ™ d 

minimising said function using co-ordinate descent optimisation; 
wherein R is a coefficient matrix of the system of linear equations; h is a vector of 
the N unknown variables; f3 is a vector containing the value of the right hand side of 
25 each equation; R(m,n) is an element of the matrix R; h(m) is the mth element of the 
matrix h; and p (n) is the nth element of the vector fi . 



The present inventors have discovered that solving a system of linear equations by 
minimising a quadratic function using co-ordinate descent optimisation offers 
30 considerable and surprising efficiency benefits. 
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According to a fourth aspect of the present invention, there is provided a computer 
processor configured to solve a system of N linear equations in N unknown 
variables, comprising: 
5 storage means for storing an estimate value for each unknown variable; 

storage means for storing coefficients of each unknown variable in each 
equation; 

storage means for storing a scalar value d\ 
initialising means for initialising each estimate value; 
10 computing means configured to process each estimate value by determining 

whether a respective predetermined condition is satisfied, and to update the estimate 
if and only if the respective predetermined condition is satisfied, said computing 
means being configured to repeatedly process each estimate value until each estimate 
value is sufficiently close to an accurate value of the respective unknown variable. 

15 

According to a fifth aspect of the present invention, there is provided a computer 
processor configured to solve a system of N linear equations in N unknown 
variables, comprising: 

storage means for storing an estimate value for each unknown variable; 
20 storage means for storing coefficients of each unknown variable in each 

equation; 

storage means for storing a scalar value d; 
initialising means for initialising each estimate value; 
computing means configured to: 
25 (a) attempt to update each estimate value using a scalar value d, 

(b) update the scalar value d if no updates are made in step (a); 

and 

(c) repeat step (a) and step (b) until each estimate value is 
sufficiently close to an accurate value of the respective unknown 

30 variable. 
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According to a sixth aspect of the present invention, there is provided a multiuser 
receiver device for obtaining data transmitted by a plurality users, the device 
comprising: 

a plurality of filters, each filter being arranged to filter out a spreading code 
5 used by a respective user; 

equation solving means to find a solution h of a system of linear equations of 
the form Rh = p where R is the cross correlation of the spreading codes used by the 
plurality of users, and fi is a vector containing the filter output signals; and 

means for obtaining the transmitted data using a solution provided by the 
1 0 equation solving means; 

wherein the equation solving means: 

(a) stores an estimate value for each value of the solution h ; 

(b) initialises each estimate value to a predetermined value; 

(c) for each estimate value: 

15 (i) determines whether a respective predetermined condition is 

satisfied; and 

(ii) updates the estimate if and only if the respective 
predetermined condition is satisfied; and 
(d) repeats step (c) until each estimate value is sufficiently close to an 
20 accurate value of the respective unknown variable. 

According to a seventh aspect of the present invention, there is provided a multiuser 
receiver device for obtaining data transmitted by a plurality of users, the device 
comprising: 

25 a plurality of filters, each filter being arranged to filter out a spreading code 

used by a respective user; 

equation solving means to find a solution h of a system of linear equations of 
the form Rh = /? where R is the cross correlation of the spreading codes used by the 
plurality of users, and J3 is a vector containing the filter output signals; and 

30 means to obtain the transmitted data using a solution provided by the equation 

solving means; 
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wherein the equation solving means: 

(a) stores an estimate value for each unknown variable; 

(b) initialises each estimate value to a predetermined value; 

(c) attempts to update each estimate value using a scalar value d\ 

(d) updates the scalar value if no updates are made in step (c); and 

(e) repeats step (c) and step (d) until each estimate value is sufficiently 



close to an accurate value of the respective unknown variable. 



According to an eighth aspect of the present invention, there is provided a method for 
10 generating filter coefficients for use in an echo cancellation apparatus, the method 
comprising: 

(a) generating a cross correlation matrix R containing the cross 
correlation of first and second signals and ; 

(b) generating an auto correlation vector containing an autocorrelation 
15 of the first signal; and 

(c) determining a vector h for which Rh = ft , said vector h containing 
the said filter coefficients; 

wherein the vector h is determined by solving the system of equations Rh = by: 

(d) storing an estimate value for each element of the vector h; 
20 (e) initialising each estimate value to a predetermined value; 



(g) repeating step (f) until each estimate value is sufficiently close to an 
accurate value of the respective unknown variable. 

According to a ninth aspect of the present invention, there is provided a method for 
30 generating filter coefficients for use in an echo cancellation apparatus, the method 
comprising: 



(f) 



for each estimate value: 



25 



(i) determining whether a respective predetermined condition 
is satisfied; and 

(ii) updating the estimate if and only if the respective 
predetermined condition is satisfied; and 
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(a) generating a cross correlation matrix R containing the cross 
correlation of first and second signals; 

(b) generating an auto correlation vector fi containing an autocorrelation 
of the first signal; and 

(c) determining a vector h for which Rh = fi , said vector h containing 
the said filter coefficients; 

wherein the vector h is determined by solving the system of equations Rh = by: 

(d) storing an estimate value for each unknown variable; 

(e) initialising each estimate value to a predetermined value; 

(f) attempting to update each estimate value using a scalar value d\ 

(g) updating the scalar value if no updates are made in step (f); and 

(h) repeating step (g) and step (h) until each estimate value is sufficiently 
close to an accurate value of the respective unknown variable. 

1 5 Brief description of drawings 

Embodiments of the present invention will now be described, by way of example, 
with reference to the accompanying drawings, in which: 

Figure 1 is a flow chart of an algorithm for solving linear equations in accordance 
20 with the present invention; 

Figure 2 is a graph showing how the value of the solution vector changes as the 
algorithm of Figure 1 solves a system of equations; 

25 Figure 3 is a graph showing how the error of the solution vector varies as the 
algorithm of Figure 1 solves the system of equations; 

Figure 4 is a graph showing how the number of passes through the system equations 
varies in dependence upon the bit number of the solution vector elements being 
30 considered in the algorithm of Figure 1 ; 
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Figure 5 is a graph showing the number of updates to the solution vector carried out 
for each bit as the algorithm of Figure 1 solves the system of equations; 



Figure 6 is a flow chart showing a variant to the algorithm of Figure 1 ; 

5 

Figure 7 is a MATLAB code fragment implementing the algorithm of Figure 6; 

Figure 8 is a flow chart showing a variant to the algorithm of Figure 6, where values 
of the solution vector are constrained between upper and lower bounds. 

10 

Figure 9 is flow chart showing how the algorithm of Figure 1 can be adapted to solve 
equations having complex coefficients and complex solutions; 

Figure 10 is a flow chart showing a method for determining the next course of action 
15 in a step of Figure 9; 

Figure 1 1 is a flow chart showing a variant to the algorithm of Figure 10; 

Figure 12 is a MATLAB code fragment implementing the algorithm of Figure 11; 

20 

Figure 13 is a schematic illustration of a device configured to solve linear equations 
in accordance with the invention; 

Figure 14 is a schematic illustration of the equation solving microprocessor of Figure 
25 13; 

Figure 15 is a schematic illustration showing block R of Figure 14 in further detail; 
Figure 16 is a schematic illustration showing block h of Figure 14 in further detail; 

30 

Figure 17 is a schematic illustration showing how the block h of Figure 16 is 
modified when equations having complex valued solutions are to be solved; 
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Figure 18 is a schematic illustration showing block Q of Figure 14 in further detail; 

i 

Figure 19 is a schematic illustration showing the minimisation block of Figure 14 in 
5 further detail; 

Figure 19a is a MATLAB code fragment showing a variant to the algorithm of 
Figure 7; 

10 Figure 20 is a schematic illustration of a CDMA receiver device in which algorithms 
of the present invention may be applied; 

Figure 21 is a schematic illustration of an echo cancellation apparatus in which the 
algorithms of the present invention may be applied; and 

15 

Figure 22 is a schematic illustration showing part of the apparatus of Figure 21 in 
further detail. 

Description of preferred embodiments 

20 

A method for a system of solving linear equations is now described. A system of 
linear equations can be expressed in the form: 

Rh = £ (1) 

25 

where: R is a coefficient matrix of the system of equations; 
h is a vector of the unknown variables; and 

(3 is a vector containing the value of the right hand side of each equations 
30 For example, the system of equations (2): 
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15x + 5^-2z = 15 
5* + 1 ly + 4z = 47 
-2x + 4y + 9z = 51 



can be expressed in the form of equation (1) where: 
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To solve the system of equations, it is necessary to find values for x, y, and z of h 
which satisfy each of the three equations. 

10 In operation, algorithm uses the matrix R and the vectors h and f5 as set out above, 
together with an auxiliary vector Q. The vector h is initialised to a predetermined 
initial value (see below) and updated as the algorithm proceeds until its elements 
represent the solution of the equations. 

15 For a system of N equations in N unknown variables, the vector h has length N and 
the matrix R is of size N x N. 

Referring to Figure 1, at step SI the vectors h and Q are initialised. The vector h is 
initialised such that all its elements are set to '0\ The vector Q is initialised to 
20 contain the negative of the equivalent position of /? . That is: 

Q = -/? (4) 

Therefore, when working with system of equations (2), Q is initialised in accordance 
25 with equation (5): 
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Q = 



-15 
-47 
-51 



(5) 



The algorithm maintains three counter variables p, m and it, a parameter N which 
represents the number of elements in the solution vector (and also the number of 
5 equations), a parameter M which represents the number of bits used to represent each 
element of the solution vector h, a parameter Nit which represents the maximum 
number of iterations through which the algorithm can pass for a particular value of 
m, a variable Flag which is used to indicate whether or not the solution vector has 
been updated, and a constant H, the purpose of which is described below. 

10 

Some of these variables are initialised at step S2 and step S3 of Figure l.p,m and it 
are all initialised to zero. N, M y and Nit are set to the values described above which 
can either be chosen by a user or hard coded into the algorithm. Selection and 
initialisation of H is described below. 

15 

Operation of the algorithm can be summarised as follows. Each bit m of all elements 
p of the solution vector h is considered in turn. As described below, for each bit an 
element of the vector Q is compared with various conditions and the result of this 
comparison determines whether or not further processing is applicable. This further 
20 processing comprises an appropriate update of the element p of the solution vector h 
corresponding to the element being considered and updates of all elements of the 
auxiliary vector Q. 



When it is determined that further processing for that element is not appropriate (for 
25 the current bit), the next element is considered. When each element has been 
considered for that particular bit, all elements of the solution vector are considered 
for the next bit in turn, and updated appropriately. This process continues until all 
elements have been considered for all bits. If the total number of iterations for any 
one bit reaches a predetermined limited the algorithm again ends. The algorithm is 
30 described in further detail below. 
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At step S3, the value of m is incremented to e l\ Thus, the algorithm is now 
considering the first bit of each element in the solution vector h. it is set to 0 to 
indicate that no iterations have yet taken place for the current value of m. 

5 

At Step S4, a step size parameter d is calculated according to the equation: 



10 where m and H are the parameters described above. 

H is a value greater than or equal to the magnitude of the maximum value which is 
expected for any value of the solution vector. That is, the algorithm considers only 
solutions lying between -H and +H. 



As will be described below, setting d in accordance with equation (6) allows each bit 
of each value of the solution vector h to be considered in turn. 

At Step S5 of Figure 1, the variable it is incremented, and the variable Flag is set to 
20 '0'. p (the current element of the solution vector under consideration) is incremented 



Having performed the necessary incrementation and initialisation, the algorithm can 
begin processing elements of the matrix and vectors, in an attempt to solve the 
25 equations. 

At step S7, the following operation is performed: 



d = 2~ m H 



(6) 



15 



at Step S6. 




(7) 



30 



where 



argmin = - 



1 , if |QQO < -Q(p) a Q(p) < -R(p, p) • -j 

2, if j- Q(p) < Q(p) a -Q(p) < -R(p, p) ■ - 



3, if j- R(p, p) • | < Q(p) a -R(p, />) • | < -Q(^) 



(8) 



The value of arg is assessed at the decision block of step S8. 

5 If arg = 1, the element of the solution vector under consideration, that is h(p) is set 
according to equation (9) at step S9. 



h(p) = h(p)+d 



(9) 



10 The auxiliary vector Q is then updated such that all values of Q are set according to 
equation (10) at step S10. 



Q(r) = Q(r) + dR(p,r), Vr : 1 < r < N 



(10) 



15 If arg = 2, the element of the solution vector under consideration, that is h(p), is set 
according to equation (1 1) at step SI 1. 



h(s>) = h(p)-d 



(11) 



20 The auxiliary vector Q is then updated such that all values of Q are set according to 
equation (12) at step S10. 



Q(r) = Q(r) - <JR(p, r),\fr : 1 < r < N 



(12) 



25 If arg = 1 or if arg = 2, Flag is set to '1' at step S13. 
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If arg = 3, no update is made to any element of the solution vector h or the auxiliary 
vector Q, and Flag is not updated. 

Having made the updates set out above, a decision block at step S14 checks the 
5 condition of equation (13); 

p=N (13) 

If p is not equal to N (i.e. all elements of the solution vector h have not yet been 
10 considered), control returns to step S6 and p is incremented. This process continues 
until all entries in the solution vector h have been considered, and h and Q are 
updated in the manner set out above. 

If p is equal to N (step SI 4), a check is made to determine the value of Flag (step 
15 S15). 

Flag is initially set to '0' at step S5, and only updated (to be equal to T) if entries of 
the solution and auxiliary vectors are amended by steps S9 and S10, or steps SI 1 and 
S12. Thus, if Flag = 1, it can be deduced that a change was made to at least one 

20 element of h (i.e. one h(p) value) and all values of Q, for the current iteration it. 
Therefore, assuming that the total number of iterations it has not exceeded the limit 
set by NzY (checked at SI 6), p is reset to 4 0' at step SI 7, control returns to step S5, 
and the current bit is again processed for each element p of the solution vector h. 
This is because further processing of each element of h for the current value of m 

25 may result in further updates. If the total number of iterations has reached the limit 
set by Nit (step SI 6), the algorithm exits. 

If it is the case that Flag = 0 (step 15), it can be deduced that no updates have been 
made to any elements of the solution vector h or the auxiliary vector Q for any value 
30 of p (that is any element of the solution vector h). Given that further iterations of 
steps S5 to S15 will result in no changes to the elements of h (given that neither d, h 
nor Q have changed), a check is made to determine whether or not all bits m have 
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been considered (step SI 8), by comparing the current value of m with the total 
number of bits M. If it is the case that m=M, i.e. all bits have been considered, there 
is no work for the algorithm to do, and the algorithm exits at step SI 9. 



5 If it is the case that all bits have not been considered, p is reset to 0 at step S20, and 
control returns to step S3, and the algorithm processes the next bit of the solution 
vector entries. 

In the preceding discussion, it has been explained that entries of the solution vector 
10 h are processed for each bit of the solution vector entries, starting with the most 
significant bit. However, it can be seen from the preceding discussion, that at all 
steps the entire value of an element of h is used for update. However bit wise 
processing is in fact achieved because following each increment of m (step S3) a new 
value of d is created at step S4. Given that each increment of m will result in the 
15 value of d being divided by two (given the presence of the expression 2~ m in the 
expression for d), and given that d is used to update both h and Q, after an update of 
d the next most significant bit is then updated. 

It has been described that the value of H represents a value greater than or equal to 
20 the magnitude of the maximum value of the solution vector elements. In setting H it 
is desirable to ensure that it is a power of two. That is, H is set according to equation 
(14). 

H = 2 q (14) 

25 

where q is any integer (i.e. positive, negative or zero) 

When //is set in this way, the expression for d set out in equation (6) becomes: 



30 a (15) 

d - 2 9 ~ m 
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Thus, when H is chosen in accordance with equation (14), the value of d can be 
updated without multiplication or division, simply by appropriate bit shift operations. 
This is particularly advantageous, because microprocessors can typically carry out bit 
shift operations far more efficiently than multiplication or division. 



The application of the algorithm described with reference to Figure 1 to the system 
of linear equations (2) set out above will now be described. 



R and /? are initialised as described above: 



R = 



15 5 -2 
5 114 
-2 4 9 



(16) 



15 
47 
51 



(17) 



h and Q are initialised in the manner described above: 



h = 



0 
0 
0 



(18) 



Q = 



-15 
-47 
-51 



(19) 



Variables are initialised as follows at step S2 and step S3: 
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m=0 
it = 0 
p = 0 



M = 8 
N = 3 
Nit = 0 



H is in this case set to 16, i.e. q = 4 in equation (14). 

m is incremented such that m=l, and it is set to c 0' (step S3). A value of d is 
5 computed according to equation (15) and in this case: 



it is incremented to *V and Flag is set to '0' at step S5. p is incremented to 1 at step 
10 S6. At step S7, the following expression is evaluated 



arg = argmin{- 15,15,-60} 

Therefore, it can be seen from equations (8) and (22) that arg =3. The decision block 
at step S8 therefore passes control to step SI 4, where the condition p=N is checked. 
15 In this case p=l and N =3, therefore p is not equal to N, and therefore control passes 
to step S6, with no changes having been made to the elements of h or Q. 

Step S6 then increments p to be equal to 2 (i.e. the second element of the solution 
vector is being considered) 




(22) 
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Step S7 computes: 



arg = arg min |q(2) ,-Q(2) -R(2,2) • || 

arg = arg min-! - 47,47,-1 !•-> 
arg = arg min{- 47,47,-44} 



(23) 



Therefore, arg=l. 

S8 the appropriate decision outcome is chosen, and at step S9, h is set as 



At step 
5 follows: 



h. 



"0" 




0 




"0" 


0 

_0_ 


becomes y 


Q + d 
0 




8 






0 



(24) 



At step S10, all values of Q are set by adding the current values to the values of the 
10 second row (since p=2) of R multiplied by d: 



(25) 



"-15" 




"-15 + t/-R(2,D" 




" -15 + 8-5 " 




" 25 " 


-47 
-51 


becomes ^ 


_47 + rf-R(2,2) 
-51 + rf-R(2,3)_ 




-47 + 8-11 




41 






-51 + 8-4 




-19 



At step 



S13, Flag is set to '1' to show that h and Q have been updated. 
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At step SU P is still not equal to N (p=2, N=3) and therefore control returns to step 
S6, where p is incremented to 3. 

At step S7 the following is computed: 



20 



20 



arg = arg min |Q(3),-Q(3),-R(3,3) ■ - 



arg = arg min<! - 1 9,1 9,-9 - — 



(26) 



arg = arg min{- 19,1 9,-36} 

Therefore arg=3. No updates are made, and the condition of step S14 is checked. In 
this instance /?=N=3, and the condition is therefore true. The value of Flag is then 
checked at step SI 5. Flag was set to T at step S13 while p was set to 2, and 
5 therefore the condition of step S15 evaluates to false. 

The condition of step S16 is then checked. Given that iY=l, and Nit=lQ, the 
condition of step S16 returns False, and control passes to step SI 7 where p is reset to 
'0' and then to step S5, where it is incremented and Flag is reset to '0'. 



p is incremented at step S6, and at step S7, the following expression is considered 



15 It can be seen from equation (27) that arg =3. The decision block therefore passes 
control to step SI 4, where the condition p=N is checked. In this case p=l and N=3 9 
therefore p is not equal to N, and control therefore passes to step S6, with no changes 
having been made to the elements of h or Q. 

20 p is incremented to 2 at step S6 and the following expression is considered at step 
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(27) 



arg = arg min {25,-25,-60} 



S7: 



,-R(2,2).f} 



arg = argmin|Q(2)-Q(2) 

arg = argmin|4 1,-41 -1 1 ~j (28) 
arg = arg min{4 1,-4 1,-44} 



It can be seen from equation (28) that arg =3. The decision block therefore passes 
control to step S14, where the condition p=N is again checked, p is still not equal to 
N 9 and therefore control passes to step S6, with no changes having been made to the 
5 elements of h or Q. 

The following expression is then considered: 



arg = argmin|Q(3) -Q(3) ,-R(3,3).|| 



arg = arg minj- 1 9,1 9,-9 • || (29) 
arg = arg min{- 19,1 9,-36} 



10 Again, arg=3, and no updates are made. In this case the condition of step S14 returns 
true, and the condition of step S15 also returns true, given that no updates where 
made during the last pass through all elements of h, and consequently Flag is still set 
to '0'. 



15 The condition of step SI 8 is then checked to determine whether all bits of the 
solution vector entries have been considered. In this case m = 1 and M= 8, therefore 
step SI 8 returns false, p is set to '0' at step S20, and control passes to step S3 where 
m is incremented to 2 and it is reset to *0\ 

20 At step S4 d is set where m is equal to 2, and therefore d=4. Note that as discussed 
above, d has been halved. At step S5 it is incremented to 1 and Flag is set to '0'. 

p is incremented to 1 at step S6. At step S7, the following expression is considered: 
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arg = argminjQ(l) -Q(l) -R(l,l)-|j 

arg = argmin|25,-25 -15 ~| 
arg = argmin{25 -25 -30} 



(30) 



Therefore arg=3, and no update takes place. The algorithm continues as described 
above, and p is incremented to 2 at step S6. At step S7 the following expression is 
5 evaluated: 



arg = argmin |Q(2) ,-Q(2) -R(2,2) • - j 

arg = arg min|4 1,-41,-1 1 -~| 
arg - arg min{41,-4 1,-22} 



(31) 



In this case arg=2. The decision block of step S8 therefore directs control to step SI 1. 
Step Sll updates h by subtracting the current value of d (d=4) from the second 
1 0 element of h, as shown in equation (32): 



"0" 




0 




"0" 


8 


becomes ^ 


8 d 




4 






0 




0 




0 



(32) 



Step S12 updates Q by subtracting the row of R (given that p=2) multiplied by d 
1 5 from the current values of Q as set out in equation (33): 



Q: 



25 
41 
-19 



becomes 



25-</-R(2,l) " 
41-rf-R(2,2) 
-19-rf-R(2,3) 





" 25-4-5 " 




5 




41-4 11 




-3 




-19-4- 4 




-35 



(33) 



Flag is set to 1 at step SI 3. 
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Another iteration is then carried out, wherein p is set to 3 at step S6, and the 
following expression is considered at step S7: 



arg = argmin|Q(3),-Q(3) -R(3,3) -|| 



arg = argmin j- 35,35 -9 ■ |J (34) 
arg = arg min{- 35,35,-18} 

Therefore, arg is set to 1, the decision block of step S8 directs control to step S9, and 
steps S9 and S10 set h and Q as set out below: 
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h: 



"0" 




0 




~o" 


4 


becomes 


4 




4 


0 




_0 + d 




_4_ 



(35) 



" 5 " 




5 + </-R(3,l) 




" 5 + 4--2 " 




"-3" 


-3 


becomes ^ 


-3 + rf-R(3,2) 




-3 + 4-4 




13 








-35 




_-35 + </-R(3,3)_ 




-35 + 4-9_ 




1 

_ _ 



(36) 



Given that p=N and Flag=l,p is reset (step SI 7) and control passes to S5. For each 
15 of the three values ofp steps S6 to S14 are executed. On each occasion, arg=3 and no 
updates take place. The individual expressions considered by step S7 during each 
iteration are not set out in full, although they can be readily deduced from the 
information presented above. 

20 Given that Flag=0, the algorithm continues for the next value of m. 



Execution of the algorithm then continues in the manner outlined above, for each bit 
of the solution vector elements in turn. 
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The values of h after each iteration through the solution vector elements are set out 
below. It should be noted that iteration number referred to here is equivalent to a 
cumulative iteration count instead of the bit by bit iteration count it described above. 



Iteration: 0 12 

1st element: 0 0 0 

2nd element: 0 8 8 

3rd element: 0 0 0 



3456789 10 

0 0 0 0 1 1 1 1 

4422222 2 

4444555 5 



Figure 2 is a graph showing the value of each solution vector element after each 
iteration. A first line 1 represents changes to the first solution vector element, h(l), a 
second line 2 represents changes to the second solution vector element h(2), and a 
third line 3 represents changes to the third solution vector element h(3). 

Solving the set of equations set out at (2) above in a conventional way yields 



h = 



(37) 



Thus, it can be seen that the algorithm effectively solves the system of equations 
after seven passes through the solution vector elements. 



The error in the values of h after each iteration is shown below: 



Iteration: 0 1 2 3456789 10 

lstelement: -1 -1 -1 -1 -1 -1 -1 0 0 0 0 

2nd element: -26 6 2 2 0 00000 

3rd element: -5 -5 -5 -1 -1 -1 -1 0 0 0 0 
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These values are plotted in the graph of Figure 3, where a first line 4 represents 
changes to the error of the first solution vector element, h(l), a second line 5 
represents changes to the error of the second solution vector element h(2), and a third 
line 6 represents changes to the error of the third solution vector element h(3). It can 
be seen that errors diminish as the algorithm proceeds. 

Figure 4 is a graph showing the number of iterations carried out for each bit, i.e. the 
value of it when processing of each bit m has been completed. Figure 5 shows the 
number of updates made to the solution vector h and auxiliary vector Q for each bit 
m. 



As described above, the auxiliary vector Q is updated as the algorithm progresses. 
The values of Q after each update of the vector Q and the solution vector h are set 
out below: 



Update 0 1 

1st element -15 25 

2nd element -47 41 

3rd element -51 -19 



2 3 4 5 6 

5 -3 -13 2 0 

-3 13 -9 -4 0 

-35 1 -7-9 0 



As a further example, consider the system of equations set out below: 
15jc + 5j/-2z=15 

5x + ll>; + 4z = -37 (38) 
-2x + 4.y + 9z = -55 



In this case the accurate solution of the equations is: 



h = 



1 

-2 
-5 



(39) 
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It should also be noted that the value H is now equal to 256. Other parameters of the 
algorithm remain unchanged. 



5 The value of the solution vector after each pass of the algorithm is set out below. 
Again, it can be seen that the algorithm correctly solves the system of equations. 



Iteration: 012345 6 7 8 9 10 11 12 

1st element: 000000 0 0 0 0 01 1 

2nd element: 00000 0 0 0 0 -2-2-2-2 

3rdelement: 0 0 0 0 0 -8 -8 -8 -6 -6 -6 -5 -5 



(40) 



10 Figure 6 shows a variant to the algorithm described above with reference to Figure 1. 
Many steps of Figure 6 are identical to steps of Figure 1. Such steps are identified by 
like reference numerals in both Figure 1 and Figure 6. Only steps which differ from 
Figure 1 are described in further detail below. 

15 In Figure 6, Step S4 of Figure 1 is replaced by step S21. It can be seen that in 
addition to setting d in accordance with equation (6): 

d = 2~ m H (6) 

20 a two element array delta is established as follows: 

delta[l]=d (41) 
delta[2] = -d (42) 

25 In Figure 1 step S8 determines the value of arg (i.e. 1, 2 or 3) and chooses a different 
action in dependence upon the value. In Figure 6, step S8 is replaced with a single 
comparison at step S22: 

arg<3 (43) 
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„ the — refoma fa.se, if can be defermined ft* arg=3 (gtven - 
amy ever take vahres of T, V, and '3-). Therefore in accordance wth fine afgonfhm 
of Figure 1. Figwe 6 shows that if ore condition is false, no updates to h or Q are 
made, and control passes to step S14 as in Figure 1. 

,f the condition of step S22 of Fignre 6 is satisfied, it can be deduced that arg - 1 or 
arg . 2 In Figure 1, if arg = 1, steps S9 and S10 update h and Q using expresstons 

, - o stem Sll and S12 update h and Q using 

including d. Similarly, if arg - 2, steps Ml ana v 

1 0 expressions including -d. 

fa to variant of the a,gorithm shown in Figure 6, the steps S9 and S, 1 of Figure 1 
are rep.aced with a single step S23, and similar* steps S10 and 812 are replaced 
wttt, a single step S24. Bout of steps S23 and S24 invoWe the array deUa and more 
15 particularly the element arg of the array delta. 

Step S23 updoes b according to equation (44) and step S24 updates Q according to 
equation (45). 

(44) 

20 h{p) = h(p) + delta(xg) 

Q(r) = Q(r ) + j e / to (arg) • R(p,r), Vr : r = 1,..., N 

Given that the first element of delta contains d and the second element contains -d it 
can be seen that equations (44) and (45) correctly correspond to equivalent 
25 operations of the algorithm of Figure 1. 

Flgure 7 shows MATLAB® source code for implementing the algorithm illustrated 
in Figure 6. 

30 Figure 8 is a flow char, showing a mrmer variant to the algorithm described above. 
Tbe afgorifhm of Figure 8 is intended ,o be used where upper and lower bounds are 
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t0 be enforced upon elements of the solution vector n. For example, it may be known 
the algorithm if only values in this range are considered. 

Figur e 8 illustrates an algorithm where all values h(p) of the solution vector h are 
constrained with the bounds of equation (46): 

(46) 

h x <h(p)<h 2 

Before execution of .he aigorithm constants », and h are estabhshed The 
establishment of these constants is not illustrated in the flowchart of F.gure a. Steps 
of Fig »re 8 which are identica! to steps of Figure 6 are identified by hire reference 

described in further detail below. 
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Having updated Up) at step S23, a test is executed a. a decision block of step S2 
ensure that new.y set value of *) lies between fire bounds specified in equafon 
(46) If the new va.ue of *) satisfies equation (46), flte algorithm proceeds «o 
update the auxihary vector Q a. S24 and to update Fla g a. step 8» as desenbeo 
20 above with reference to Figure 6. 

If step S25 determines that the newly se, value of h(p) does not satisfy equation (46) 
h(p) is again updated a. step S26. This updating reverses tire update of step 823. such 
2 «p) is equal to the value before execution of step S23. Wore has no 
25 been updated (because the at.enrp.ed update did no, satisfy equation (46)), F.ag ,s no. 
set to «1% and execution of the algorithm continues at step S14. 

to yet a further embodiment of the present invention, the constants t, and fc of the 
embodiment of Figure 8 axe replaced by vectors h, and * bom 
30 length equal to the solution vector b. The vector h, contains a lower bound for each 
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element of the solution vector, and the vector h 2 contains an upper bound for each 
element of the solution vector. Thus, the condition of step S25 becomes: 



h,0>)<h(/>)<h 2 Q>) 



(46) 



The embodiments of the present invention descnbed above are concerned with the 
application of the algorithm to systems of equations which have real valued 
solution, The present invention is also applicable to the solution of systems of 
equations having complex valued solutions, and the application of the invention to 
such systems of equations is described below. 

Consider the system of equations set out below: 
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(a, + a 2 i)x + (b, + b 2 i)y + (c, + c 2 i)z = A,+ A 2 i 
(</, + d 2 i)x + (e, + e 2 i)y + (/, + = + ^ 
(g, + g 2 *)x + (K + h 2 i)y + + j 2 Q* = C, + C 2 i 



(47) 



where each of the 
follows: 



unknown variables x, y and z is a complex number defined as 



20 and 



x = x { + x 2 i 
y = yi + y 2 i 
z = z.+ z 2 i 



(48) 



From 



equation (1) above, the system of equations (47) can be expressed as follows: 
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R = 



fl ,+a 2 j b x +b 2 i c x +c 2 i 
d x +d 2 i e x +e 2 i A+fi 1 
gx+gi* hi* 11 !* Ji+Jii 



(49) 
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h = 



jc, +x 2 i 



z, + z 2 i 



(50) 



A x + A 2 i 
B x +B 2 i 
C, + C 2 i 



(51) 
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To solve the system of equations, a matrix A, and vectors b, and c are created from 
the data set out above. A is a 2N by 2N real-valued coefficient matrix, b is a real- 
valued solution vector of length 2N, and c is a real- valued right hand side vector of 
length 2N, where N is the number of unknown variables (i.e. A/=3 in this example). 

A, b, and c are set as described below: 



A(2w - l,2/i - 1) = Re{R(m, «)} 
A(2m,2n) = Re{R(w,n)} 
A(2m - l,2n) = Im{R(/w,n)} 
A(2w,2« - 1) = - Im {R(m, «)} 



(52) 
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b(2«-l) = Re{h(«)} 
b(2«) = Im{h(n)} 



(53) 



c(2«-l) = Re{>9(n)} 
c(2«) = Im{/?(«)} 



(54) 



Where Re{} is a function returning the real coefficient of a complex number, and 
20 Im{} is a function returning the imaginary coefficient of a complex number. 



Thus, in the case of equations (47), A, is set as follows: 
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A = 



«1 


°2 


by 


b 2 




c 2 




«1 


-b 2 


bi 






d x 


d 2 


*x 


e 2 


/. 


/l 


-d 2 


dx 


-*j 


*x 




fx 


gi 


gi 


K 


h 2 


» h 


1 


-Si 


g\ 


-K 


*. 


-h 





(55) 



The present invention is often used to solve normal equations. Wherenormal 
equations are solved, the matrix A is such that: 



b x =d, 
b 2 =~d 2 

fx-K 
f 2 =' h 2 

C x =gx 
C 2 = -g 2 



a, >0 
e, >0 

a 2 =0 
e 2 =0 
i, =0 



(55a) 



The vectors b and c are set as follows: 
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„ avi „ g crea*, <ne vecrora a ano b ano *. ma** A se, out above, .he — * 
sowing real va,ued eauauons se, out above W i<h reference ,0 Frgures ., 6 ana 8 can 



be used where: 
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R = A (58) 
h=b 



= c 



VaH.es of the solution vector h se, by the agoriUtm can then be usee to define 
Mb ,be real and imaginary components of.be complex numbers *, y and a, and 
5 create the complex valued solutions. 

* alternative embodiment of the present invention, the a.gorithms described ab.ve 
JZm - ,at the algorithm operates direct, on R, h and as set out tn 
anions (49, to (5.) above. Th- modifications are now described with reference 
10 to the flow chart of Figure 9. 

The algorithm used to solve conations involving comp.ex numbers is b*ed upon that 
of Fi ie ., bu, steps S7 to S,2 are replaced by the steps shown m Ftgure 9. S*p 
S27 oTFigure 9 replaces step S7 of Figure 1, step S2S of Figure 9 replaces «p » 0 
15 rl . 1 steps S29 to S36 replace steps S9 to S» of Figure ,. Steps 

SHarl'idenbca, , me erpvfien, steps of Figure , and are shown in dotted hues. 

They are included in Figure 9 for the sake of clarity. 

Referring to Figure 9, at step S27, the following expression is evaluated: 
20 j-y 

arg .a J gmm{Re(QW>,-Re(Q(p)},Im(Q(rt>,-Im(Q(p)},-R(P.P) I j < 5 » 



where 
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argmin = 



1 , i fj mi »(Re{Q(p)},-^(Q(P)!.'»^»-- ta ' Q(rt! '- R(P ' P) I ) = Re(Q(rt) l 
2,if(min(R={Q(p)},-Re ( Q(p)), ta (Q(P)>.-^»(<^(P».- R (''•''>•f ) = - Re(Q(P), 

3 , i f{ m i„(Re { Q(p)},-Re(Q(p)), ta (Q(rti.-^«Wl-- R( ''-'' ) ! ) = Im!Q(P)} l 

4,if{min(Re(Q(p)},-Re{Q(rt).ta(Q(P)).- ]m «^> ) '- R( ' , -' ,) f ) = - Im{Q(P)> 
5,if mta(Re{Q(p)},-Re{Q(P)i.In.(Q(P)).-ta(Q(P)!.-R(i'.rt l) = - R( ^ rt 1 



(60) 



and Re(i and Im{) are as defined above. 



The value of arg set at step S27 is checked a, the decision block of S28. This 
amines the update (if any) to be tnadeto the elements of the solution vector h and 
the auxiliary vector Q. 

10 If arg is 1, .hen the element p of the solution vector h currently under consideration is 
updated as follows at step S29: 



h(p)=h(p)+d 



(61) 



Equation (61) means that the real part of h(p) is increased by d. 



Ml elements of the auxiliary vector Q are also updated, in accordance with equation 
15 (62) at step S30. 

Q{r ) = Q(r) + dR(p,r)yrA<r<N 



If arg is 2, then the element p 



of the solution vector h currently under consideration is 
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updated as follows at step S31: 
h(p)=h(p)-d 



(63) 
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Equation (64) means that the real part of h(p) is decreased by d. 

All elements of the auxiliary vector Q are also updated, in accordance with equation 
(62) at step S32. 

Q(r) = Q(r) - dR(p, r),Vr :l<r<N (64) 

5 

If arg is 3, then the element p of the solution vector h currently under consideration is 
updated as follows at step S3 3: 

h(p)= h(p)+vd (65) 

Equation (65) means that the imaginary part of h(p) is increased by d. 

10 All elements of the auxiliary vector Q are also updated, in accordance with equation 
(66) at step S34. 

Q(r) = Q(r) + d •/ - R(p 9 r), Vr : 1 < r < N (66) 

If arg is 4, then the element p of the solution vector h currently under consideration is 
15 updated as follows at step S3 5: 

h(p)= h(p)-id (65a) 

Equation (65a) means that the imaginary part of h(p) is decreased by d. 

All elements of the auxiliary vector Q are also updated at step S36, in accordance 
with equation (66a). 

20 Q(r) = Q(r) - d • i • R(/>, r), Vr : 1 < r < N (66a) 

If arg = 1, 2, 3, or 4, Flag is set to '1* at step SI 3, control then passes to step S14 and 
the algorithm proceeds as described above with reference to Figure 1 . 

If arg = 5, no updates are made to h or Q and control passes to step SI 4, and the 
algorithm proceeds as described above. 
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Step S27 described above, requires the minimum of five values to be identified so as 
to determine the next course of action: 

min{Re{Q(p)},^ 

a method for finding this minimum, and efficiently identifying the correct action (at 
5 step S28 of Figure 9), is now described with reference to the flowchart of Figure 10 

The algorithm takes as input two values a and b. a is input at a step SI 01 and is set to 
be equal to the real part of Q(p), and b is input at a step SI 02 and is set to be equal to 
the imaginary part of Q(p). 

A decision block SI 03 determines whether or not a is greater than 0. If a is positive, 
10 a is set to be equal to the negative of its current value at a step SI 04, and a variable 
Ja is set to T at step SI 05. If a is not positive, the variable Ja is set to '0' at step 
SI 06. 

b is processed in a similar manner. Step SI 07 checks whether b is positive. If b is 
positive, it is set to the negative of its current value at step SI 08, and a variable Jb is 
15 set to T at step S109. If b is not positive, Jb is set to '0* at step SI 10. 

Thus after execution of steps S 101 to SI 10, Ja is set to T if the input value of a was 
positive, and '0' if the input value of a was not positive. Jb is similarly set. Given the 
action of steps SI 04 and SI 08, it can be seen that both a and b will now not be 
positive. 

20 At step Sill a check is made to determine whether a is greater than b. If this 
condition is true a variable Jl is set to be equal to Jb, a variable J2 is set to '0' and a 
variable c is set to be equal to b. This is accomplished by step SI 12. If the condition 
of step Sill is false, a variable Jl is set to be equal to Ja, a variable J2 is set to '0' 
and c is set to be equal to a. This is accomplished by step SI 13. 

25 At step SI 14, a check is made to determined whether or not c is greater than - 
R(p,p)d/2. If this condition returns true, a variable J3 is set to be equal to T at step 
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SI 15. If this condition is false, the variable J3 is set to be equal to '0' at the step 
S116. 

The method of Figure 10 is such that after execution of all steps, the variables Jl, J2 
and J3 are set as follows: 

J3 = 1 => no update necessary 
J3 = 0 update necessary 

J2 - 1=> update to imaginary part necessary 
J2 = 0 => update to real part necessary 

Jl = 1 => update should comprise subtraction 
Jl = 0 => update should comprise addition 

It will be appreciated that the modifications made to the algorithm of Figure 1, as 
illustrated in Figures 6 and 8 can similarly made to the algorithm of Figure 9. For 
example, the algorithm of Figure 9 can be implemented using an array delta, similar 
to that used in Figure 6, where: 

delta[l] = d; 
delta[2] = -d; 

delta[3] = d-i; (6?) 
delta[4] = -d-i; 

Having established delta in this way, the algorithm can proceed by simply adding the 
appropriate element of delta to the appropriate element or elements of h or Q as 
described above. Figure 1 1 shows this variant of the algorithm of Figure 9 where 
steps S28 to S36 are replaced by steps S37 to S39. It will be appreciated that 
additionally the array delta must be initialised and updated after each update of d in 
accordance with equation (67). This is not shown in Figure 11. Figure 12 shows a 
MATLAB program implementing the algorithm illustrated in Figure 11. 
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The de.crip.ion presented above has illustrated varions implementahons of 
a>gori.hms in aceordanee with the invention. It wi.l be appreciated that vanons 
amendments can be made to these algorithms without departing from the invention. 

For example, in the description set out above execution ends when the number of 
iterations * for any particular bit m reaches a predetermined limit MB. It wtll be 
appreciated that execution need not end in this circumstance. Instead, a tuner t may 
be set to 'IT each time m is updated, and execution can end if this timer exceeds a 
predetermined time threshold. 

The vectors h and Q need not necessarily be initialised as indicated above. Indeed, 
the initial value for h should usually be substantially centrally positioned m the range 
-H to H in which solutions are being sought, so as to obtain quick convergence. In 
some embodiments of the present invention the anxiliar, vector need not be created. 
Instead the vector p is used directly, and is updated in a manner similar to the 
manner described above for vector Q, although it will be appreciated that different 
updates will be required. Suitable updates for such embodiments of the inventton are 
set out in the derivation of the algorithm presented below. 

It has been described above that d is updated by dividing the previous value of d by 
two This is preferred because considerable benefits are achieved because 
compotations involving multiplication of d can be carried out using efficient btt sh.fi 
operations in place of relatively inefficient multiplication operations. However, .. 
will be appreciated that alternative methods for updating d may be used m some 
embodiments of ft. invention. For example, if d is updated by division by a power of 
two such as four or eight, computations can still be efficiently implemented by 
carrying out two or three bit shifts instead of me single hi. shifts required when d .s 
updated by division by two. 

, in preferred embodiments of me present invention, each value of me solotion vector 
h is represented by a fixed point binary word. This is particularly beuefic.al gtven 
that mathematical operations car, typically be carried ou. more efficiently usmg 
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fixed poin, arithmetic. Furthermore, a fixed point representation is likely to be 
aeeeptable beeause the different unknown variables are like.y to have an 
approximately equal magnitude. 

to eirenmstanees where a fixed point representation is inappropriate for the solution 
veetor values, a eonventional floating point representation ean be used. Although the 
algorithm will operate more s.owly with floating poin. values than with fixed potnt 
values, the algorithm still offers very favourable performanee as eompared wtth other 
methods for solving linear equations. 

It will be appreciated that the algorithm described above can be implemented in 
software or hardware. A software implementation will typically comprise appropnate 
source code executing on an appropriate microprocessor. For example, as shown m 
Figures 7 and 12, code can be created in MATLAB which is compiled to create 
object code which is specified in the instruction set of a microprocessor. Although 
MATLAB provides a convenient implementation for the algorithms of the invention, 
it will be appreciated that the algorithms could instead be coded in any one of the 
large number of widely available computer programming languages. 

A hardware implementation of the algorithm can either be provided by configuration 
of appropriate ^configurable hardware, such as an application specific integrated 
circuit (ASIC), field programmable gate array (FPGA), or a configurable computer 
(CC) or alternatively by a bespoke microprocessor built to implement the algorithm. 
Figures 13 to 18 illustrate a device and microprocessor configured to solve linear 
25 equations. 

Figure 13 sehematieally illustrates the general arehiteeture of a deviee configured to 
solve linear equations. The deviee oomprises a host proeessor 7 which is a general 
purpose mieroproeessor of known form configured to control operation of the dev.ce 
30 by running appropriate program code. This program code can be stored on a non 
volatile storage deviee 8 (e.g ROM or FLASH memory) and read by me genera, 
purpose microprocessor 7 when it is to be executed. Alternatively the program code 
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can be copied to a volatile memory 9 (e.g RAM) prior ,o execution. In order to solve 
teear equations nsing the method of the invention, the device comprises an equatton 
solving microprocessor 10, operation of which is controlled by the host processor 7. 
The aforementioned components of the device are connected together by a host bus 
11. 

Figure 14 illustrates the architecture of the equation so.ving microprocessor 10 in 
tether detail. The equation solving microprocessor 10 comprises a controller 12, the 
function of which is to control operation of the equation solving microprocessor 10. 
The equation solving microprocessor 10 further comprises a h-b.ock 13, which stores 
and updates the sohttion vector h of the algorithm, a Q-block 14 which stores and 
upda.es the auxiliary vector Q of the algorithm, and a R-b,ock 15 which stores the 
coefften. matrix R used by the algorithm. The equation solving microprocessor 10 
also comprises a minimisation block !6 which finds the minimum of a number of 
vatues. Components of the equation solving microprocessor 10 are connected 
together by an internal bos 17, along which control instructions and data can pass. 

Figure 15 shows the structure of R-block 15 in further detaif. The R-block 15 
comprises a storage element 18 which stores .he elements of .he coefficient matnx 
R The R-block also comprises muttiplexers 19, 20 and 21 which generate a colunrn 
address signal 22, and a row address signal 23 for reading dafa from and writing date 
to .he storage element 18 in me manner described below. Finally, me R-block 
comprises a bit-shift e.emen, 24 which can left-shift a value presented a. its input by 
a value determined by the current value of d. 

The multiplexers 19, 20, 21 generate address signals as follows. The row multiplexer 
20 and the column multiplexer 21 each have selection inpufs connected to a common 
input line 25 which con carry an initialisation signal Mt which is typically recerved 
tan the controller 12 (Figure 14). If an initialisation signal is recetved, tins 
j indicates drat date is to be written to the storage element 18 representing the 
coefficients of .he equations which are to be solved. In mis case me row mulhplexer 
20 should generate a row address 23 which is equal to the input row address 26, and 
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similar* the column multiplexer 2! should generate a eolumn address 22 wraeh ,s 
equa. to the input eo.nmn address 27. During initialisation, the row address 26 and 
cohunn address 27 wil. typically eount through e.entents of the matrix R wntmg daU 
provided ,o the storage element 18 on input line 28 to me appropriate e.emen, of the 
5 storage element 18. 

When initialisation has been completed, the InU signal wi.l no longer be received by 
me selection input of the multiplexers 20, 21, and therefore the row address ,s 
determined by the input 29 to the row multiplexer 20 and Ore eolumn address ,s 
,0 determined by .he input 30 to the column multiplexer 21. The input 30 is me value, 
of the algorithm as described above. The inpu, 29 is connected to the update 
multiplexer 19 whose output is dependent upon whether values of R are being rea^ 
for purposes of ana.ysis <e.g step S7 of Figure 1) or update of Q (e.g. steps S.O and 
S12 of Figure 1). The update multiplexer 19 has a selection inpu. 31 which mdrca.es 
,5 either update or analysis. This can conveniently be achieved, for example, by a s mgle 
bit input where '0' indicates update and T indicates analysis. 

If R is being read for purposes of analysis (i.e. the selection input 3 1 is se, 
then the element R(p,,) is required (see step S7 of Figure 1). Therefore, update 
20 multipiexer output 29 is se. to be equal to the input 32 which is equal to „, and the 
row multiplexer 20 subsequently generates a row address 23 equal top. 

If R is being read for purposes of update of the auxiliary vector Q (i.e. selection 
inpu. is se, to '0') Ffer) is require* (see steps S10 and S12 of Fignre 1). Therefore 
25 the update multiplexer output 29 is equal to the input 33 which is set to be equal to ,. 
and the row multiplexer 20 subsequently generates a row address 23 equal to r. 

In either case, the column multiplexer 21 generates a column address 22 which is 
equal to p. 

30 

From the preceding description, it can be seen that the multiplexers 19, 20 and 21, 
act together to ensure that the cc.ee, row address 23 and column address 22 are sen, 
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to the storage element 18. Having read the appropriate element of R from the storage 
element 18, the output needs to be multiplied by d regardless of whether it is being 
used for update or analysis (see steps S7, S10 and S12 of Figure 1). As described 
above, this can be achieved by a bit shift (assuming that d is a power of 2), the 
number of bits to shift being determined by the expression: 

, a (68) 
log 2 rf 

which is passed to an input 34 of the bit shift block 24. The output 35 from the bit 
shift block 24 is then correctly set as follows: 

■ c , (68a) 
dR(p,r), lfupdate 

, (68b) 
d-R(p,p), if analysis 



Figure 16 illustrates the h-block 13 of Figure 14 in further detail. The h-block 
comprises a storage element 36 for storing elements of the solution vector h, an 
update/reading multiplexer 37, and an adder-subtractor 38. 

To initialise elements of the solution vector, an initialisation signal 39 is input to the 
storage element 36. Upon receipt of this signal, all elements of the solution vector h 
are initialised to '0'. 

The update/reading multiplexer 37 sends an appropriate address 39a to the storage 
element 36. The address 39a sent to the storage element 36 is determined by the 
signal provided to a selection input 40 of the multiplexer 37. If the selection input 40 
is set to update, the address p provided at input 41 becomes the address 39a. 

If the selection input 40 is set to read (i.e. the values of h are being read to obtain the 
solution of the equations), a read address 42 input to the multiplxer 37 becomes the 
address 39. In the case of a read operation, the input read address will count through 



42 



all elements of h in turn, and each element will be transmitted to the internal bus 17 
of the equation solving microprocessor 10 in turn. 

In the case of an update operation, a single value of p is provided to the multiplexer 
5 37, and an input signal I\ is provided to the adder/subtractor 38. The signal I\ 
indicates whether the element of the solution vector h(p) is to be updated by adding d 
or subtracting d. Upon receipt of the address 39, the storage element 36 outputs the 
appropriate element to the adder/subtractor 38 along a line 44. By analysing the 
signal I\ the adder/subtractor 38 can determine whether its output should be set to 
10 h(p)+d or h(p) - d. The appropriate expression is calculated, and written back to the 
storage element 36 along a line 45. 

Figure 17 shows how the h-block 16 can be implemented when the algorithm is used 
to solve equations having complex values (for example using the algorithm of 
15 Figures 10, 1 1 or 12). It can be seen that in addition to the components illustrated in 
Figure 16, the update/read multiplexer 37 has an additional input 46 which carries a 
signal This signal is sent to the storage element 36 along a line 47 along with the 
address p when the selection input is set to update. The signal I 2 indicates whether 
the update should be made to the real part or the imaginary part of h(p). 

20 

Figure 18 illustrates the Q-block 14 of Figure 14 in further detail. It can be seen that 
the Q-block comprises a storage element 48 for storing the vector Q, an address 
multiplexer 49, a data multiplexer 50 and an adder/subtractor 51. 

25 The address multiplexer 49 has a selection input 52 which selects update, analysis, or 
initialisation mode. The address multiplexer 49 also has three data inputs, a first data 
input 53 carries the value r used in the algorithm, a second data input 54 carries the 
value p used in the algorithm, and a third data input 55 carries initialisation address 
data. 

30 

When the selection input 52 of the address multiplexer 49 is set to initialise, an 
output 56 of the address multiplexer carries the initialisation address supplied at the 
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« ,„pu, 55 of the adore, nruuip.exer. When the se.ect.on n»* « - - » 
update the output 56 earries rite va,ue , supplied a. the firs, tnpu, 53 of m address 
JJ. ^ when the se.ecrion inpu, « is se, ,o ana,ysis, * 
address mu.rip.exer 49 is se. to .he value p suppfied a. * second tnpu,* th 
address multiplexer 49. The output 56 of <he address mu.Up.exer , passed .o 
storage elemen. 48 as an address for reading or writing dala. 

Da ta is written to Are storage Cement 48 o„.y in the update and hn,ia.isa,ion modes^ 
T 1 ore .he da, ntn.fip.exer 50, has a se.ecrion input 57 having two recogntsed 

X 'e common with the se.ecrion input 52 of dte address ntu.rip.exer 49, a,.ho gh 
Zbehaviour of the data ntu.rip.exer is no, we,, defined when the se.ecrion tnpu, ts 
set to analysis). 

15 The data ntuUip.exer 50 comprises a firs, da,a input 58 carrying huriahsa rion data 
(Centers of the vector „ and a second da. input 59 carrying update data. «en 

data inpu, 58 is provided a, dte outpu, 60 of the data tnu,«p,exer 50. When rite 
HhtpuJisse, to update, data fiom the second data inpu, 59 , provtdcd a, 

20 the output 60. 

„ operation, when rite se.ecrion inpu* of both ntuUip.exers 49 50 rue se, to 
LL a sequence of addresses (counting riuough a„ Cements of me vector Q, are 

25 1 ritese addresses, the appropriate da,a is provided a, ,he first data tnpu 58 oHhe 
data mu,rip,exer 50. The da,a and addresses are provided to me 
and me storage elemen, is men se, ,o contain me appropnate Cements of the vector 



Q 



30 When rite se.ecrion inputs of ho,h muhtp.exers 49, 50 are se, to update, rite -*-r 
I provided to the storage Cement 48 hy the ntuhiplexer 49. The 
Z *e storage elemen, 48, Q(r, ts provided at an output 6, of the storage Cement, 
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and passed ,0 a first input 62 of the addcr/s„b«rac«or «• ^ 
addlbtractor is set to ,*,,) (P-iaed by the R-block > - « 
adder,sub«rac.or also includes an input 64 ca^ying a stgna. /, whob md, cates 
fite adder/subractor 5, should perfomr an addfiion or subtract™ operas 
T l i , of appropriate data a, its inputs an output 65 of the adde—or ts set 

I Given that tbe selection input 57 is se, .0 update, fire data provided at the m u, 5 

pLed to fire output 60 of the data multiplexer, and fit,. , there ,0 e s,ora 8 e 
JLm 48 where i. is written to the element of Q specified by the address 56. 

' Where the algorithm is being implemented to so.vc complex valued equations «he 
abactor comprises a further input 66 carting an signa, H whtch^ 
whedter the upda* should be appfied to the tea, par. or the tmagmary pat. of 
appropriate element of the vector Q. 

5 h ana.ysis .node, ute Q-b,ock is required to provide a curren. va.ue of QM with no 
X Una bemg needed. The data mul.iplexer 50 ,s merest — 
1 fire analysis operation. The address tnu.fiplexer 49 provrdes fire address , 

20 Q W is read front fire s.orage element 48 and provided a. ft. ou,pu. 61 of fite storage 
element 48. 

Figure .9 provides an implementation for the minimisation block 16 of Figure 14 
X m P — ton corresponds to ma. Ulustrated in Figure 10, and ts there o. 

The minimisation block comprises firs, and second converter blocks 67, 68, firs, and 
second comparators 69, 70 and first and second multiplexers 7 1, 72. 

fcmut data is provided to the firs, and second converter blocks 67, 68. The lb* 
3„ ri b 1 67 takes as input the real par. of a value «* and me seco ; d 
letter block 6S takes as tnpu. me imaginary part of the vame Q(p). The — 
,„cks 67 and68p ro videtwooutpu,s,afirs,outpu.67a, 68a indicates me stgnofme 
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input data, and second output 67b, 68b indicates the modulus of the input data. The 
modulus outputs 67b, 68b are input to the first comparator 69, which generates a 
signal I 2 which at an output 73. The signal I 2 is defined as being '0' if the value 
provided at output 67b is greater than that provided at 68b, and T if the value 
5 provided at output 68b is greater than or equal to that provided at output 67b. 

The values provided at outpus 67b and 68b are also provided to the first multiplexer 
71, along with the output value 73 created by the first comparator 69. The output 73 
of the first comparator acts as a selection input to the first multiplexer 71. The first 
10 multiplexer 71 then acts to provide the larger of the values 67b, 68b at its output 74. 

The output 74 of the multiplexer 71 is input to the second comparator 70, along with 
the value d R(p f p) at an input 75, which is computed by the R-block 16 of Figure 15. 
The second comparator 70 then produces an output 76 which is a signal I3 which is 
15 set such that I 3 is set to '0' if the input 74 greater than the input 75, and T if the 
imputer 75 is greater than or equal to the input 74. 

The second multiplexer 72 is a bit multiplexer taking as input the signs of the input 
data generated by the converter blocks 67, 68. In dependence upon the value of the 
20 signal I 2 output from the first comparator 69, the multiplexer generates an output 77 
which is a signal Ij. 

The signals //, I 2 and I 3 can together be used to determine what update (if any) should 
be made to the elements of Q and h. I 3 indicates whether or not the current iteration 
25 is successful. If I 3 is equal to '0', the element h(p) is updated, and all elements of Q 
are updated. If I 3 is equal to ' 1', no updates are necessary. 

h indicates whether the real or the imaginary part of the appropriate elements should 
be updated. If I 2 is equal to '0' the real part of the appropriate values is updated, 
30 while if I 2 is equal to ' V the imaginary part of the appropriate values is updated. // 
indicates whether the update should comprise addition or subtraction. If // is set to 
c 0', addition is used, and if // is set to ' 1\ subtraction is used. 
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Having described the structure and function of the individual components of the 
equation solving microprocessor 10 (Figure 14), operation of the processor is now 
described. 

5 

To perform initialisation, the equation solving microprocessor receives a signal to 
cause initialisation, for example from the host processor (Figure 13). The controller 
12 of the equation solving microprocessor then generates an appropriate internal 
initialisation signal which is communicated to all blocks of the equation solving 

10 microprocessor, via the internal bus 17. Upon receipt of this signal the h-block 13, 
the Q-block 14 and the R-block 15 perform initialisation as described above. The 
data required for initialisation may be located in a predefined read only memory, or 
an address where the data is to be found may be provided to the controller 12 for 
onward transmission to the appropriate block. Parameters used by the algorithm (e.g. 

15 N, M and N/Y) can either be provided to the controller, or alternatively specified 
within the microprocessor 10. 

When initialisation is complete, the controller begins executing an algorithm in 
accordance with the invention using the blocks of the microprocessor 10 to update 

20 values of h and Q as appropriate. In some embodiments of the invention data that is 
required to be passed between blocks of the microprocessor is passed directly 
between blocks, under the control of the sending and/or receiving block. In 
alternative embodiments of the invention all data is passed between blocks as 
directed by the controller 12. When the controller determines that the equations 

25 have been solved, such that the vector h contains the solution of the equations, the 
vector h can then be copied from the storage element 36 (Figures 16 and 17) to an 
appropriate location within the device, where the solutions can be used as required. 

It will be appreciated that although a specific hardware implementation of algorithms 
30 of the invention has been described above, numerous modifications could be made to 
the implementation while remaining within the scope and spirit of the invention. 
Such modifications will be readily apparent to those of ordinary skill in the art. 
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The preceding description has described algorithms in accordance with the invention, 
and hardware suitable for implementing such algorithms. The mathematical basis of 
algorithms in accordance with the invention is now described. 



To aid understanding of the invention, the known co-ordinate descent optimisation 
method for minimising a function of many variables is presented. The co-ordinate 
descent optimisation method seeks to find: 



where J is a function of many variables; and 

h is an N-dimensional vector containing the variables. 

1 5 Thus we want to find the value of h when the function J is a minimum. 

In many circumstances the elements of h have a maximum amplitude H. Therefore, 
the minimisation problem is considered for a subset of values of h defined by 
equations (70) and (71). 



5 



10 



min{J(h)} 



(69) 



20 



heU 



(70) 



where: 



25 




(71) 



where: 



h(m) are elements of the vector h and H is a known number such that H>0. 
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Define a vector e, where the ith coordinate is T, and all other coordinates are '0'. 
Let h 0 be an initial value of the solution vector h and let a 0 be an initial value of the 
step size parameter (that is d in the algorithms described above). 



5 hk is defined to be the value of the solution vector h after some number of iterations 
k, where k is a positive integer. a k is the step size parameter after k iterations. 

A vector p* is also defined, according to the equation: 
10 p k = e ik (72) 

where: 

i k =k-N(div(K,N)) + \ (73) 

15 where div (A,B) is a function whose result is defined as the integer part of result of 
dividing A by B. 

It can be seen from equation (73) that i k will repeatedly count from 1 to N as k is 
increased. This results in values of p cycling through ei to eN: 

20 

The value of the function J is calculated at the point h = h k + a k j>k. and two 
conditions are checked: 

25 

h k +a k p k e\J (75) 
J(h k +a k p k )<J(h k ) (76) 
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If the conditions of equations (75) and (76) are satisfied, then the solution vector h is 
updated as follows: 

(77) 

5 

Given that that update involves a scalar multiple of the vector p k and given also that 
p k is a vector containing only a single non-zero element, it can be seen that equation 
(77) means that h* + i will be identical to h* save for a single element. 

10 The step size parameter is not updated at this stage, as indicated by equation (78) 

(78) 

If the conditions of equations of (75) and (76) are not satisfied, the value of the 
15 function J is calculated at the point h = h* - « tP , and two conditions are again 
checked: 

J(h k -a k p k )<J(h k ) (80) 

20 

If the conditions of equations (79) and (80) are satisfied, then the solution vector h is 
updated as follows: 

(81) 

h*+i = h*- a k pic. 

Again, given that that update involves a scalar multiple of the vector p k and given 
also that p k is a vector containing only a single non-zero element, it can be seen that 
equation (81) means that h k+x will be identical to h* save for a single element. 

30 The step size parameter is not updated at this stage, as indicated by equation (82) 



25 
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(82) 

The ft* iteration is considered to be successful if either the conditions of equations 
(75) and (76) or equations (79) and (80) are satisfied. If neither of these conditions 
are satisfied, the iteration is considered to be unsuccessful. 



a k is updated according to equation (83): 



10 a k +\,= 



[Aa t ,ifi t =WAh t =h t _„ +I ( 8 3) 
I a k , otherwise 



where X is a parameter of the algorithm, and is such that 0 < X < 1 

The condition of equation (83) means that if the last N iterations involve no 
15 successful updates (i.e. the value of k has not changed), the step size parameter 1S 
updated by multiplication by X . If however there is at least one update during the 
previous N iterations, the step size parameter is not updated. It should be recalled that 
N is defined to be the number of elements in the solution vector h, and it can 
therefore be seen that h is updated only when all elements have been processed and 
20 no update has occurred. 

The method described above is generally known as co-ordinate descent optimisation. 
It is known that for some arbitrary function J, the method converges to find the value 
of h for which the function is a minimum, providing that the function J is convex and 
25 differentiable on U, providing that a 0 > 0, and 0<X<\ and providing that h 0 e U . 
This result is shown, for example, in Vasiliev, F.P, "Numerical methods for 
solutions of optimisation problems", Nauka, Moscow 1988 (published in Russian), 
the contents of which is herein incorporated by reference. 
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It is often necessary to solve the linear least squares problem. The linear least squares 
problem is concerned with the minimisation of the function J specified in equation 
(84): 

5 J(h) = |Zh-d| 2 = (Zh-d) r (Zh-d) < 84) 

with respect to an unknown vector h, where Z is a known M x N matrix, d is a 
known vector of length N, and " r denotes the transpose of a matrix. 

10 This is discussed in Sayed, A.H., and Kailath, K, "Recursive least-squares adaptive 
filters", The Digital Signal Processing Handbook, CRC Press, IEEE Press 1998, 
pages 21.1 to 21.37, which discussion is herein incorporated by reference. 

It can be shown that the minimisation of the function of equation (84) is equivalent 
1 5 to minimisation of a quadratic function. 



Equation (84) can be rewritten as: 

/(h) = h r Z r Zh - h r Z r d - d r Zh + d'd 



20 



by multiplying out the bracketed expressions of equation (84). 



Given that the purpose of the method is to minimise J with respect to h it can be 
concluded that the term d r d of equation (85) will not effect the minimisation process, 
25 and therefore can be removed from the expression without affecting the minimum 
value. Therefore minimisation of equation (85) is equivalent to minimisation of 
equation (86): 



j(h)-h r Z r Zh-h r zWZh < 86) 
A matrix R is defined according to equation (87): 

T ( 87 ) 

R = Z r Z 



30 
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A matrix p is defined according to equation (88): 



t (88) 

5 

Equation (86) can then be rewritten using R and fi as shown in equation (89): 

JW-h T Bh-h T fi-P T * (89) 



10 Simplifying this expression yields: 

J(h) = h T Kh-2h T /3 (90) 

The expression h r Rh can be rewritten in terms of the elements of h, and R as 
15 follows, using the definitions of matrix multiplication, and the matrix transpose 



operation: 

h r Rh = f j f j R(m,n)h(m)h(n) 



(91) 



Similarly the expression h T /3 can be written in form of equation (92): 

20 

h r /? = f>(n)h(nj 



(92) 



Substituting equations (91) and (92) into equation (90) yields: 



25 J(h) = R(w. n)h{m)Hn) - 2^ P(n)h(n) 

m=l «=1 B=1 



(93) 



It can be seen that equation (93) is a quadratic function of h. Thus it can be seen that 
solving the linear least squares problem is equivalent to minimisation of the function 
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of equation (93). Furthermore, it is also known that solving a system of linear 
equations of the form of equation (1): 



(1) 
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is equivalent to minimisation of the function of equation (93), for any set of normal 
linear equations. This is explained in, for example, Moon, Todd K., and Stirling, 
Wynn C: "Mathematical methods and algorithms for signal processing", Prentice 
Hall, 2000, section 3.4. "Matrix representations of least-squares problems", pages 
10 138-139. This explanation is incorporated herein by reference. 

Given that many sets of linear equations occurring the electronics and physics are 
normal linear equations, minimisation of the function of equation (93) has wide 
applicability in solving linear equations. 



The explanation presented above has set out a method for minimising a function J 
using the co-ordinate descent optimisation method. The material presented above has 
also set out the relationship between the minimisation of equation (93) and a set of 
linear equations of the form of equation (1). 



The present inventors have surprisingly discovered that applying the known co- 
ordinate descent method to the minimisation of equation (93) provides a particularly 
efficient method for solving linear equations. 

25 Minimisation of the function of equation (94) is considered: 
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(94) 




This minimisation process finds values for the elements of h which minimise the 
function J(h). The matrix R and the vector ft are known. It is known that the 
function of equation (94) is convex and differentiable. This is shown in Vasiliev, 
F.P.: "Numerical methods for solutions of optimisation problems", Nauka, Moscow 
5 1988 (published in Russian), page 345, which explanation is incorporated herein by 
reference. Therefore, as explained above, the co-ordinate descent optimisation 
method can be used to find the minimum value of the function J, 

During operation of the co-ordinate descent optimisation method, the following 
10 expressions are computed: 

U(h k ) = J(h k +a k e ik )-J(h k ) (95) 

and 

15 AJ(h k ) = J(h k -a k e ik )-J(h k ) (96) 

It can be seen that equation (95) relates to the condition of equation (76) set out 
above, while equation (96) relates to the condition of equation (80) set out above. 

20 

The inequality of equation (97) is also considered: 

A/(h,)<0 (97) 
25 It can be recalled that the values of the vector e, has elements defined by: 



fl ifi = i* 

KL=i (98) 
* [0 otherwise 

Also, it is known that the matrix R is symmetric, given that the system of equations 
30 is normal. Therefore, substituting equation (94) into equation (95) yields: 
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(99) 



where h (k) (m) are elements of the vector ta k . 
5 If a vector Q is defined as: 

Q (l ) = -2^(i) + 2Xh w (m)R(m,0 (100) 

m=l 

then equation (99) can be written as: 

10 

A/(h k ) = ^[Q(0 + « i R(M)] (101) 

Given that a k > 0 , from equation (101), equation (95) can be rewritten as: 

(102) 

a*R(M)<-Q$ 

15 

Similarly, equation (96) can be rewritten as: 

(103) 

a*R(u)>Q$ 

20 An auxiliary vector Q k is defined as the vector Q after the ft* iteration. If the (^l)th 
iteration is not successful, then elements of h and Q can be updated as follows: 

(104) 

hft+i = hk 

25 and 

(105) 

Qa+i = Qk 
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If the (fc+l)th iteration 
follows: 



is successful, then elements of h and Q are updated as 



( 106 > 

5 h^\i)=h (k) (i)±a k 

q(*+D ( n ) = <f>(n) ± 2a k R(n,i),n = 1,..., # (107) 

10 The vector h can be initialised to a vector h 0 where 

(108) 

h 0 (n) = 0,n = \,...,N 

15 

Then, from equation (100), elements of Q are initialised as follows: 

(1° 9 ) 

Qo(n) = -2P(/i;,n = l,-.^ 

20 . . f 

That is , each clement of Q is set to be the negative of the eorrespondtng element of 

p multiplied by two. 

Thus, the solution veotor h is initialised to contain all V values, while the auxiliao- 
25 vector Q is initialised to be the negative of the vector /J . 

As described above, multiplication operations may be avoided by setting H in 
accordance with equation (110): 

(110) 

30 H = 2"^ 

where M b is a positive integer, and P is any integer. 

a 0 is initialised to be -y and X is initialised to be i . 

The multiplications described above can then be replaced by bit shift operation, 
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The algorithms described thus far have used an auxiliary vector Q which is intialised 
in accordance with equation (4): 

(4) 

Q = -fi 

However, some embodiments of the invention use fi itself as an auxiliary vector. In 
such embodiments of the invention, the auxiliary vector update rule of equation (107) 
above becomes: 

P ^(n) = P (k) (n) +a k R(.n,i),n = \,...,N ( nl > 
Similarly, the inequalities of equations (102) and (103) become: 

15 a*R(u)<2£(z) (U2) 

and: 

(113) 
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a*R(M)> 

Figure 19a shows MATLAB source code for an algorithm which uses fi in place of 
Q. This algorithm is based upon that of Figure 7, but the code has been amended as 
set out above. 

The step size parameter a, is updated if N consecutive iterations are not successful. 
At every update of a k , a k is decreased by a factor of two. The algorithm described 
above is therefore referred to as the dichotomous coordinate descent algorithm for 
the solution of linear equations. 

30 From the description set out above, it can be observed that the algorithms of the 
invention solve linear equations by minimisation of an appropriate quadratic 
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function. It has been explained above that it is known that such minimisation can be 
employed to solve normal linear equations. However, it should be noted that the 
present invention is not limited simply to normal linear equations, but is instead 
applicable to a wider class of linear equations. 

In the methods described above, the elements of the solution vector h are analysed in 
a predetermined order (i.e. from element T to element AO- However, it will be 
appreciated that elements of the solution vector h can be analysed in any convenient 
manner. For example, the values of h can be sorted on the basis of some 
corresponding auxiliary value (e.g. a corresponding element of the vector Q), and 
elements of the solution vector h can then be processed in that order. For some 
applications, ordering elements of the vector h in this way will provide a more ra P1 d 
convergence, although this must of course be balanced against the computational 
cost of sorting the elements of h. 

It has been explained above that the present invention can be usefully applied in any 
application in which it is necessary to solve linear equations. Two such applications 
are now described. 

In a multiuser Code Division Multiple Access (CDMA) communications system, a 
plurality of users transmit data using a common collection of frequencies. A narrow 
band data signal which a user is to transmit is multiplied by a relatively broad band 
spreading code. Data is then transmitted using this broad band of frequencies. Each 
user is allocated a unique spreading code. 

A receiver needs to be able to receive data transmitted by a plurality of users 
simultaneously, each user using his/her respective spreading code. The receiver 
therefore needs to implement functions which allow the spreading code to be 
removed from the received data to yield the originally transmitted data. Typically 
filters are used to extract the spreading code to obtain the transmitted data. It should 
be noted that the process is complicated by interfering signals from multiple users, 
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and also from different propagation paths which may be followed by different 
signals. 

Figure 20 illustrates a receiver 80 suitable for use in a CDMA communications 
5 system. The receiver comprises a receiver circuit 81 including an antenna (not 
shown), an analog to digital converter 82, a bank of filters 83a, 83b and 83c, an 
equation solving circuit 84 and a decision circuit 85. 

Spread spectrum signals are received by the receiver circuit 81, and converted into 
10 digital data by the analog to digital converter 82. Digital data is then passed to all of 
the filters 83a, 83b, 83c. Each filter of the bank of filters 83a, 83b, 83c relates to a 
unique user, and has filter coefficients selected to match the spreading code of that 
user. It will therefore be appreciated that in practical embodiments a receiver will 
include more than three filters as is illustrated in Figure 20. 

15 

If a single signal is transmitted at any one time, and this signal travels between a 
sender and the receiver 80 by a direct path, the output of the filters alone should 
provide the data which the sender intended to transmit. However, in the more likely 
and more complicated situation where interference between signals occurs, the 
20 outputs of the filters alone will be inconclusive. However, it is known that solving an 
equation : 

( 114 ) 

Rh = j0 
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where R is the cross correlation matrix of the spreading sequences of all users; 
p is a vector containing the filter output values; and 
h is a solution vector 

will allow the originally transmitted data to be obtained.. 
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In general, for a system involving N users, there will be N filters, and the vector 
P will therefore have length N, and the matrix R will have size NxN. 



The cross correlation matrix R can be defined as 



R = S T S 



(115) 



where the matrix S contains the spreading codes, and is defined as follows: 



S = 



5,(1) 5,(1) 

5,(2) s 2 (2) 



(116) 



s x {m) s 2 (m) ••• s k (m) 



where 5j(z) denotes the ith element of the spreading code for user j\ 

As has been described above linear equations of the form shown in equation (114) 
can be solved using an algorithm in accordance with the invention. Therefore, the 
invention provides a novel multi user receiver apparatus, in which the equation (114) 
is solved as described above, thereby achieving the considerable performance 
benefits provided by solving equations in accordance with the invention. 

The equation solver 84 of Figure 20 can either be implemented by means of a 
computer program carrying out the method of the invention executing on an 
appropriate microprocessor, or alternatively by means of hardware, such as that 
described above with reference with Figures 14 to 19. 

The equation solver provides a vector h as output, and this is input to the decision 
circuit 85, which then determines the nature of the transmitted data. 
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It will be appreciated that the cross correlation matrix described with reference to 
equations (115) and (116) is merely exemplarily. Cross correlation matricies can be 
created in a variety of different ways which will be known to one skilled m the art. 
Regardless of how the cross correlation matrix is formed, a system of equations 
5 (1 14) is created which can be solved using methods in accordance with the present 
invention. 

It will also be appreciated that in addition to the components illustrated in Figure 20 
and described above, a CDMA receiver may require other components to function 
10 properly. The selection and use of these components will be readily apparent to those 
of ordinary skill in the art. 

The algorithms of the invention can also be employed in adaptive filtering 
applications such as echo cancellation in a hands free communications system. 

A system of interest in illustrated in Figure 21. A signal travels along input line 86 
and is output through a loudspeaker 87. A further signal such a human voice (not 
shown) is passed to an input of a microphone 88. It is desirable that the signal at the 
output 89 of the microphone 88 contains only the human voice signal, and none of 
20 the signal output by the loudspeaker 87. However, in practice, some of the signal 
output by the loudspeaker 87 will be received by the microphone 88 such that the 
microphone output 89 comprises a combination of the human voice signal and part of 
the loudspeaker output signal (referred to as "the echo signal"). It is desirable to 
remove the echo signal present in the microphone output 89. 

25 

As shown in Figure 21, the echo cancellation apparatus comprises a filter 90, which 
is configured to provide an estimate 91 of the echo signal. This estimate 91 is 
subtracted from the microphone output signal 89 by a subtracter 92 . Therefore, if the 
echo is accurately estimated, an output 93 of the subtractor will be equal to the 
30 human voice signal. 
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The echo cancellation apparatus comprises a filter coefficient setting circuit 94, 
which takes as inputs a signal 86 which is input to the loudspeaker and the signal 89 
which is output from the microphone. The circuit 94 should generate coefficients to 
allow the filter 90 to accurately model the echo signal. 

5 

Figure 22 shows the filter coefficient setting circuit 94 in further detail. It can be seen 
that the circuit comprises an auto correlation element 95, a cross correlation element 
965 and an equation solver 97. The auto correlation element 95 finds the auto 
correlation of the signal 86. The cross correlation element 96 finds the cross 
10 correlation between the signal 86 and the signal 89. It is known that when a auto 
correlation matrix R and a cross correlation vector /?have been generated, the 
optimal filter coefficients h can be found by solving the system of equations: 

Rh = /J (117) 

15 

where the auto correlation matrix R is generated according to the equation: 

T-\n-m\ 

R(m, n) = Yj + \ n ~ H)> ^^n:\<m<N A\<n<N (118) 

t=\ 

where t=l, ...,rare discrete time moments; 

20 

and the cross correlation vector is generated according to the equation: 

25 where x is the loudspeaker input signal 86 and y is the microphone output signal 89. 

An echo cancellation system operating in the manner described above is described in 
US5062102 (Taguchi). 
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Having generated a system of equations of the form of equation (117), algorithms in 
accordance with the present invention can be used to solve linear equations to 
determine a solution vector h containing optimal filter coefficients. Therefore, 
referring back to Figure 22, the equation solver 97 can either be a microprocessor 
executing code in accordance with one of the algorithms described above, or 
alternatively hardware suitably configured to implement a suitable algorithm. 

It will be appreciated that although this application of the algorithm has been 
described with reference to echo cancellation, it is widely applicable in all cases 
where an adaptive filter is required, and where solving a system of linear equations 
yields appropriate filter coefficients. A suitable example system in which the 
invention could be beneficially employed is described in WO 00/38319 (Heping). 

Applications of the invention to CDMA receivers, and echo cancellers have been 
described above. However, it will be appreciated that many other applications exist 
which can benefit by the improved efficiency with which linear equations can be 
solved in accordance with the invention. For example, the invention can be used in 
tomographic imaging systems, where a large system of linear equations is solved to 
generate an image. 

Although the present invention has been described above with reference to various 
preferred embodiments, it will be apparent to a skilled person that modifications lie 
within the scope and spirit of the present invention. 
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